CONDITION ESTIMATES FOR PSEUDO-ARCLENGTH CONTINUATION * 
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Abstract. 

We bound the condition number of the Jacobian in pseudo arclength continuation problems, and we 
quantify the effect of this condition number on the Hnear system solution in a Newton GMRES solve. 

In pseudo arclength continuation one repeatedly solves systems of nonlinear equations F{u{s), A(s)) = 
for a real- valued function u and a real parameter A, given different values of the arclength s. It is known that 
the Jacobian of F with respect to a; = {u, A) is nonsingular, if the path contains only regular points and 
simple fold singularities. We introduce a new characterization of simple folds in terms of the singular value 
decomposition, and we use it to derive a new bound for the norm of F~^. We also show that the convergence 
rate of GMRES in a Newton step for F{u{s), A(s)) = is essentially the same as that of the original problem 
G{u, A) — 0. In particular we prove that the bounds on the degrees of the minimal polynomials of the 
Jacobians F^ and Gu differ by at most 2. We illustrate the effectiveness of our bounds with an example from 
radiative transfer theory. 

Key words. Pseudo- Arclength Continuation, singularity, GMRES, singular vectors, eigenvalues, rank- 
one update 
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1. Introduction. Numerical continuation is the process of solving systems of nonlinear 
equations G{u, A) = for various values of a real parameter A. Here u : i?^ — i? is a 
real- valued function and G : —* . An obvious approach for implementing numerical 

continuation, called parameter continuation [6,9,15], traces out a solution path by repeatedly 
incrementing A until the desired value of A is reached. In each such iteration, the current 
solution u is used as an initial iterate for the next value of A. Although parameter continu- 
ation is simple and intuitive, it fails at points {u, A) where the Jacobian Gu is singular. In 
this paper we consider singularities which are simple folds. 

The standard way to remedy the failure of parameter continuation at simple folds is 
to reparameterize the problem by introducing the arclength parameter, s, so that both u 
and A depend on s. This idea, known as pseudo- arclength continuation [6,9,15], implements 
parameter continuation on F{u{s), X{s)) = with s as the parameter instead of solving 
G{u, A) = with A as the parameter. Thus pseudo-arclength continuation requires that the 
Jacobian of F be nonsingular. It is known that F^ is nonsingular at simple folds and 
points where Gu is nonsingular [9]. 

Our first goal (§ EI) is to quantify this nonsingularity. To this end we provide a new 
characterization of simple folds in terms of the singular value decomposition (SVD) of Gu- 
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From the SVD, we derive a new bound for ||-Far^||2- This bound can be used to hmit the 
arclength step in Newton's method. As a byproduct we obtain a refinement of Weyl's 
monotonicity theorem [19] for the smallest eigenvalue of a symmetric positive semi-definite 
matrix f ^3.H) . 

We also examine (21) how the conditioning of affects the convergence of the inner 
GMRES [20] iteration in a Newton-GMRES solver [1,2,11,12]. We show that the eigenvalue 
clustering of the Jacobian in the reformulated problem is not much different from that 
of the Jacobian G„ in the original problem. To be precise, the upper bound on the degree 
of the minimal polynomial of F^ differs from that of Gu by at most two. This implies [3, 14] 
that the convergence of GMRES as a linear solver in a Newton step does not slow down 
when parameter continuation is replaced by pseudo-arclength continuation. 

At last we illustrate our findings with a numerical example from radiative transfer 
theory. These numerical results, combined with our bounds, support the use of pseudo- 
arclength continuation in solution paths that contain simple folds. 

2. Background. We briefiy review theory and algorithms for solving numerical con- 
tinuation problems G{u,\) = 0, where X E R, u : ^ R and G : R^+^ R^ . We 
discuss parameter continuation ^2.11 and pseudo-arclength continuation in ^2.21 We use the 
abbreviations 

^dG ^dG 

OU OA 

2.1. Simple parameter continuation. Parameter continuation [6,9,15] is the simplest 
method for solving G{u, A) = 0. The idea is to start at a point A = Xinu and solve G{u, A) 
for u, say by Newton's method. Use the solution uq as the initial iterate to solve the next 
problem G{u, A+ dX) = 0. Algorithm paramc below is a simple implementation of parameter 
continuation from Xinu to Xend = Xinu + ndX where n denotes the maximum number of 
continuation iterations. 



paramc(M, G, Xinu, Xend, dX) 
Set A = Xinit, Uo = u 
while A < Xend do 

Solve G{u, A) = with Uq as the initial iterate to obtain Ui 

Uo = Ml 

X = X + dX 
end while 



While parameter continuation appears to be a reasonable method for solving G{u, A) = 0, 
it fails at points that violate the assumptions of the implicit function theorem. Such points 
of failure are called singular points. 

Definition 2.1. A singular point is a solution (mo,Ao) to G{u,X) = for which 
Gu{uo, Ao) is singular. 
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In order to understand why parameter continuation fails at singular points, we recall 
the implicit function theorem [9, 18]. The norm || ■ || denotes the Euclidean norm, and 
C^{Q) denotes the space of k times continuously differentiable functions from an open subset 

n c to i?^. 

Theorem 2.2. implicit Function Theorem:/ Let be an open subset of R^^^ and 
let G G C'^(n) for some integer k > 0. Let Gu and G\ be Lipschitz continuous in Cl, the 
closure offl. If 

• (mo, Ao) e n, 

• G(mo, Ao) = 0, 

• Guiuo, Ao) is nonsingular, 

then there are p > and e > such that there is a unique solution 

V G Bp{uo) = {u I ||n — noil < p} 

of G{u, A) = for all A G (Aq — e, Aq + e) and u G i3p(uo). Furthermore, v is a k times 
continuously differentiable function of X. 

If the assumptions of the implicit function theorem are satisfied then Newton's method 
converges q-quadratically, as shown below. We will use the following definition of quadratic 
convergence for Newton's method. 

Definition 2.3. Let C be a sequence and let x* G R^ . We say that x„ x* 
q-quadratically as n ^ oo, if Xn x* and if there is K > such that 

II * 1 1 ^ iy II * 1 1 2 

I 7T,— |— 1_ j JV il^ j • 



The following corollary presents conditions under which Newton's method applied to 
G{u, A) = converges q-quadratically. 

Corollary 2.4. Let the assumptions of Theorem \2.^A hold. Then there is 5 > 
which depends only on ||G~^(mo, Ao)|| and the Lipschitz constants of Gu and u, such that if 
|A — Ao| < 5 then Newton's method with initial iterate uq converges q-quadratically to the 
solution u* . 

Proof. Define the Lipschitz constants 
||m(A) - u{p)\\ < 7„|A - p\, ||G„(m, A) - Guiv,p)\\ < 7g(||m - f III + |A -/u|). 
According to [11] 

II _ *ll 1 

^I'<27g||G.-H%,Ao)|| 

so choosing 

1 

^27.7g||G-iK,Ao)|| 

completes the proof. □ 
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Theorem 12.21 and Corollary 12.41 suggest that points (m, A) for which Gu is singular may 
cause the loss of uniqueness in the solution to G{u, A) = as well as the failure of Newton's 
method, and therefore failure of Algorithm paramc. Thus we need a continuation method 
that does not fail at singular points. 



2.2. Pseudo-arclength Continuation. Pseudo-arclength continuation [6,9,15] avoids 
the problems of Algorithm paramc at singular points by using an arclength parameterization. 
The curve in Figure for instance, has a singularity with respect to the parameter A. If we 
choose arclength s as the parameter A, and x = [u^, X)'^ in place of u, we can compute the 
curve with simple parameter continuation. The curve in Figure l^!T] has a simple fold, which 
is the singularity of interest for this paper. Formally, a simple fold is defined as follows. 

Definition 2.5. A solution (mq, Aq) of G{u,X) =0 is a simple fold if 

• dim{Ker{Gu{uo, Xq))) = 1 and 

• Gx{uo,Xo) ^ Range{Gu{uo,Xo)). 

To develop a pseudo-arclength continuation method, we assume that x depends smoothly 
on s. Then one can differentiate G{u, A) = with respect to s and obtain 

(2.1) ''^"'W-^W' ^G„..G.A^O. 

ds 

Equivalently, one can differentiate G{x) = and obtain G^x = 0. Here, x denotes the 
derivative with respect to s. Because the norm is the Euclidean norm and s is arclength, 

(2.2) ||xf = ||^(||2 + |j,|2 = 1. 

Since we introduced a new parameter s, we must add an equation to G{u, A) = so that 
the number of equations equals the number of unknowns and we have a chance of obtaining 
a nonsingular Jacobian for the reformulated problem. Hence, we work with the extended 
equations 

(-) -(-)^(4:i))-(o)- 

The normalization equation A/" = is an approximation of ()2.2|) where 
(2.4) J^{x,s) = x^{x - xo) - {s - sq) = 0. 



This equation says that the new point on the path lies on the tangent vector through the 
current point Xq. 

Given a known point (xo,so), the pseudo-arclength continuation method increments 
arclength hj ds = s — Sq, and solves ()2.3p with the normalization ()2.4|) by Newton's method 
with initial iterate Xq. Algorithm psarc is a simple implementation of pseudo-arclength 
continuation. 
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pSarc(M, F, Send, d\) 

Set s = 0, xq = {uq, Xq)^ 
while s < Send do 
Approximate x 

Solve F{x, s) = with xq as the initial iterate obtain xi 

Xq = Xi 

s = s + ds 
end while 



Since pseudo-arclength continuation is just simple parameter continuation applied to 
F with s as the parameter, Corollary 12.41 gives conditions for the convergence of Newton's 
method in pseudo-arclength continuation. 

Corollary 2.6. Let the assumptions of Theorem \2.iA hold for F . Then there is 5 
which depends only on \\F~^{xq, Sq)\\ and the Lipschitz constants of F^ and x, such that if 
\s — Sol < S then Newton's method with initial iterate xq converges q- quadratically to the 
solution. 

The proof of Corollary 12.41 shows that the step in arclength is bounded by 



2^,^F\\F~'ixo,So)r 

where and 'jp are Lipschitz constants for x and F, respectively. Therefore a bound on 
IIF^T^II is an important factor in bounding the arclength step. In the next section we present 
the main result of this paper, a new bound on ||F^7^||. 

3. Nonsingularity of F^. For a solution xq = (uq, Aq) to G{u,X) = 0, we present an 
upper bound on ||F^^(a;o, So)|| in the case that 

• ^^(^0, Xq) is nonsingular or 

• {uq, Xq) is a simple fold of G{u, X) = 0. 

In order to derive the bound, we introduce a new characterization of simple fold, which is 
based on the singular value decomposition of Gu. We prove the bound in ^3.21 In ^3.11 
we refine Weyl's monotonicity theorem for the smallest eigenvalue of a symmetric positive 
semi-definite matrix, which we need for the proof. 
Let 

Gu{u,X) = UEV^ 
be a singular value decomposition (SVD) of Guiu, X) where 

S = diag{ai, (72, ... , ctjv), cti > a2 > ■ ■ ■ > cr^ and Un = Ucm 

where cat is the last column of the N x N identity matrix. The trailing column mjv of f/ is a 
left singular vector associated with the smallest singular value ajsf. Since the singular values 
are continuous functions of the elements in Gu{u, A), they are also continuous in A. If 

O'N-l > (J > 
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for all {u, A) then the nullity of Gu{u, A) is at most one. If in addition = then un spans 
the left nullspace of Gu{u, A). From the direct sum 



Ker{Gl) © Range{Gu) = 

we see that Gx{uo, Aq) is not in the Range{Gu) if and only if G\um 7^ 0. Hence we have a 
new, equivalent definition of simple fold. 

Definition 3.1 (Simple Fold via SVD). Let (mo,Ao) he a solution ofG{u,X) = 0, and 
let un be a left singular vector of Gu{uo, Aq) associated with a^. 

Then Aq) is a simple fold if 

• dim{Ker{Gu{uo, Xq))) = 1 and 

• u%Gx{uo, Ao) ^ 0. 

Continuity of G^un implies that there is a > such that for all {u, A) 

/ 2 I ,2 gap \ . ^ ^ 
max ajy, \u,yGx\ —j > a > 0, 

V gap + r/ 

where 

gap = cr^_i - a^, and ^ = \u%Gx\ + - UNu1f)Gx\\. 

Theorem 3.2. Let fl be the closure of an open subset Q G i?^^^, and let G be con- 
tinuously differentiable in fl. Let xq = Aq) in Q be a solution to G{uo,Xo) = 0, and 
M{xq, So) = with \\xq\\ = 1. Let r > be such that \\GuUo + GxXo\\ < r. 

Assume that for all {u, A) in Cl there exists a > such that 



(jAr_i > 0, max < alf, (ulfGxf — ^-7-77 \ > « 
I gap + 



where 

gap = o-^_i -crfi, ^ = \unGx\ + ||(/ - UNuJf)Gx\\. 

Ifr < a, then for allx = {u, A) in Cl, the smallest singular value cr^i^{Fx) of the Jacobian 
Fx of F{x, s) is bounded from below with 



We postpone the proof of Theorem 13.21 until ^3.11 in order to derive an auxiliary result 
first. 

3.1. Lower Bound for the Smallest Eigenvalue. We derive a lower bound for the 
smallest eigenvalue of the rank-one update A + yy^, where A is a real symmetric positive 
semi-definite matrix of order A^, and y is a real x 1 vector. 

Let /5i > . . . > (3n be the eigenvalues of A. Weyl's monotonicity theorem [19, Theorem 
(10.3.1)] implies bounds for the smallest eigenvalue of A + yy^'- 



pN < X^iniA + yy^) <(3n-i. 
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Intuitively one would expect that Amin(A + yy^) is larger if y is close to an eigenvector of 
(3iq. We confirm this by deriving lower bounds for Amin(A + yy"^) that incorporate the angle 
between y and the eigenspace oi [3^. 

Theorem 3.3. Let A he an N ^ N real symmetric positive semi-definite matrix, un 
an eigenvector of A associated with /Sn, \\un\\ = 1; O'lT'd y a real N x 1 vector. Set 
Vn = UNy- Then 

(3.1) Amin(^ + yy'^) > max \ I3n, yl, 



gap + ^2 



where gap = (3n-i - Pn and ^ = \yN\ + ^ lbP - Iz/ivP- 
Proof. We first show that 

(3.2) + yy^) > min{/3^ + y^^^, (3n-i%} 

gap + 4 4 

is lower bound for Amin(^ + yy^) = niin||a;||=i x'^{A + yy^)x. 
Let 

/A \ 
A = u\ ■.. 

be an eigendecomposition of A, and x be any real vector with ||x|| = 1. Partition 

so that ^ = \yN\ + \\y\\- Then 

x^{A + yy^)x > f3N-i\\xf + (3nxI + {y'^xf. 

If > \yN\/i then 

x^{A + yy^)x>{(3N-iyl)/e, 

which proves the second part of the bound in ()3.2|) . 

If ||x|| < \yN\li then \yjq\ — \\x\\^ > 0, and it makes sense to use \xiy\ > 1 — in 

\y'^x\ = \yNXN + f^x\ > lyNXwl - \\y\\ > lywl - \\xU- 

Hence 

x^{A + yy^)x > + (3^x1 + {y^xf > + + (gap + ^^)\\xf - 2^\\x\\ {y^l- 

This is a function of ||x|| which has a minimum at ||x|| = |yAr|4^/(gap + 4^^). Hence 

x^(A + yy^)x >(3n + y^— 

gap + ^2 

which proves the first part of the bound in ()3.2|1 . 
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With the help of ()3.2|) we now show the desired bound ()3.ip . Weyl's theorem [19, 
Theorem (10.3.1)] imphes Amin(A + yy^) > Pn, which proves the first part of the bound in 
(|3.ip . For the second part of the bound in (|3.1|) . we use the fact that the eigenvalues of A 
are non-negative, hence Pn-i > gap and 

pN-i ^ gap 



^2 gap + ^2 

Substituting this into (j3.2p gives the second part of the bound in (j3.ip 

mm(^ + yy )> mm{fjN + yjy yN—p^\ 

gap + r r 

^ . , 2 gap 2 gap . 2 gap 

gap + gap + e gap + 

□ 

The quantity gap in Theorem 13.31 is the absolute gap between the smallest and next 
smallest eigenvalues. The theorem shows that Amin(^ + yy^) is likely to be larger if y has a 
substantial contribution in the eigenspace of (3^- The bound in Theorem 13.31 is tight when 
y is a multiple of un- That is, if Im^tI/I = \\y\\ then Amin(^ + yy^) = minj/^Tv + WyW^, Pn-i}- 

Now we are in a position to complete the proof of Theorem 13.21 

3.2. Proof of Theorem 13.21 Define the residual r = GuUq + GxXq and form 

TP jpT _ ( Gu G\ \ / '^o \ _ f C^uG^ + G\G\ r 
" ' ^--^ \J\Gl Ao, 



- - Any VGT Any "V 1 



The eigenvalues of F^F^ are the squares of the singular values of F^. Applying Theorem! 
to GuGl + GxGl with A = GuG^, y = Gx, Pn = cr^r, Pn-i = cr%-i and gap = a%_^ - a% 
shows XminiGuG'^ + GxGJ^) > a. Hence we can write 

GuGl + GxGl ^V' p^pT ^i^E, 



1, 

where < r max l|. If r < min{a;, 1} then < 1, / + E' is nonsingular, and 

4. Newton-GMRES and Eigenvalue Clustering. This section discusses the perfor- 
mance of the inner GMRES iteration in the context of continuation with a Newton-GMRES 
nonlinear solver. Theorem 13.21 gives conditions under which the Jacobian matrix F^. of the 
reformulated problem is uniformly nonsingular. This implies GMRES is a practical candi- 
date for making the linear solve in Newton's method when implementing pseudo-arclength 
continuation. While the results in the previous section address conditioning, they do not 
directly translate into the performance of iterative methods [7,11,21], especially in the non- 
normal case. However, we can go further to see that the eigenvalue clustering properties of 
the matrix F^ do not stray far from those of Gu- 
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Suppose the eigenvalues of Gu are nicely clustered (in the sense of [3,14]). Even in the 
singular case, this would mean that the zero eigenvalue of Gu is an "outlier". We seek to 
show that adding the row and column does not significantly increase the number of outliers, 
and that we can then use the estimates in [3, 14]. 

One approach is to use the paradigm of [13]. The idea is that 

(4.1) Gu = I + K{u) + E 

where is a low-rank operator, say of rank p, and E is small. We then want to write 
in the same way, and then compare the number of outliers by comparing the ranks of the 
i^'-terms. 

Assume that E is small enough so that the eigenvalues oi I — K are "outliers" in the 
sense of [3]. Since the degree of the minimal polynomial of / — is at most p + 1, we have 
a bound for the sequence of residuals {r;} of the GMRES iteration of the form 

(4.2) ||rp+,.|| <C||E||iro|| 

where p < p + 1 GMRES iterations are needed to kill the contribution of the outlying 
eigenvalues. 

Theorem 14. II states that that the spectral properties of F^ are similar to those of Gu- 
Theorem 4.1. Let the assumptions of Theorem ]!^. Up hold. Assume that (|4.H) holds with 
rank{K{u)) = p. Then there is K,{u) having rank at most p + 2 such that 

\\F, - I - IC{u)\\ < \\E\\. 



Proof. We write [13] 

Fx = I{N+i)x{N+i) + (^^T ) + ( 

The range of 




IS 



Range{K) 






+ span < > + span 




and hence the rank of /C is at most p + 2. U 

So, while the eigenvalues may change, we have not increased the degree of the minimal 
polynomial of the main term {K vs /C) beyond p + 3. Hence, the methods of [3] can be 
applied to obtain a bound hke (j4.2|) with p < p + 3. 
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5. Example: Chandrasekhar H-Equation. We now present an example of a solution 
path containing a simple fold. The equation of interest is called the Chandrasekhar in- 
equation [4,11,17] from radiative transfer theory: 



Hiu) ^""^ 



(5.1) H{f,) = l + ^H{f,)[\ ^ , 

I Jo + u 

The goal is to compute the norm of the solution to Equation ()5.1|) for various natural 
parameter values c. That is, we compute 

\\H\\i= f\(u,c)du 
J 

as a function of c. Integrating 1)5.11) with respect to fi yields 



\H\\ =1 I j^ H{ii)H{ii)iidiidv ^ ^\\H\? 

2ioio /i + i^ 4 ^' 



and so 

(5.2) \m.^'^^ 



c/2 

Equation ()5.2|) tells us two interesting things. First, there can be no real solutions of the 
if-equation for c > 1, so there must be a singularity at c = 1, or else the implicit function 
theorem would tell us that we could continue past c = 1. Secondly, the ± gives us a hint 
that there may be two solutions, at least for < c < 1 (and there are!). 

Figure |0] is a plot of ||i^||i against c. Notice how the curve bends around when c = 1, 
and how there are two solutions for each < c < 1. In fact, we are witnessing a simple fold 
at c = 1. 



Fig. 5.1. \\H\\i as a function of c 
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5.1. Simple Fold at c = 1. For the if-equation, it is possible to compute the singularity 
analytically. Write the iJ-equation as 



G{H,c){^,)=H{f,)-ll-y 



c /■i/iif(z/) (iz/\ ^ 



IJ,+ u 

Taking the Frechet derivative of F in the direction of w yields 



Gh{H, c)w{fi) = w{fi) - ^ 2 = '^W - 7;H{fi) / 



1-f 



^fiH{iy)diy\ 2 J o fx + u 

/i + z/ J 



Let c = 1, then ()5.2|) implies that 

f^H(fi) dfi = 2 
J 



and therefore 



Hence if = fiH{fi), 



Gh{H, 1)0 = 0, 



and we have shown directly that Gh is singular at c = 1. 

One can apply Perron- Frobenius theory [8, 16] to show that the null space of Gh has 
dimension one, and hence is spanned by (j). The singularity at c = 1 is a simple fold because 

G^{H, 1)(/.) = = H\^^){H-\^^) - 1) 

2 /i + z/ 

is not in the range of Gh- To see this note that 



G,{H,l){fi) = H\fi){H-\fi) -1) <0 

and vanishes only at /i = 0. The null space of Gjj is the span of H^^, which is strictly 
positive. Hence Gc is not orthogonal to the null space of Gjj. 

One can also show that Gh is nonsingular for all c 7^ 1 by an argument even more 
tedious than the one above [10]. 



12 DICKSON et al 

5.2. Smallest Singular Values. As a demonstration of the result in § El we calculate 
the smallest singular value of the Jacobian matrix associated with the augmented system for 
the if-equation with each continuation iteration. In the language of §121 we find <Jmin{F(^H,c)) 

for various c where F(^h,c) denotes the Jacobian of ^ J^{^h'c\) ^ ' ^^^^^^ shows that the 

smallest singular value of -F(//,c) for each c stays away from zero keeping -F(_ff,c) nonsingular, 
even at the simple fold (c = 1). The pseudo-arclength code used here uses a direct LU 
factorization of the Jacobian for the linear solve in Newton's method. The step in arclength 
is fixed at (is = .5, and we use a a secant predictor [9]. The integral is discretized with the 
composite midpoint rule and 200 nodes. The singular values are calculated using Matlab's 
svd command. 



Fig. 5.2. o'niin(^(ff.c)) OS a function of c 




).1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0. 



5.3. Computation with //-equation. The consequences of the remarks in § 0] are 
that for a problem like the //-equation, which is a nonlinear compact fixed point problem, 
the number of GMRES iterations per Newton step should be bounded. One must take this 
expectation with a grain of salt because as one moves along the path, the norm of the solution 
increases, and so the number of outliers may increase slowly. The observations we present 
illustrate this. 

We use a Newton-GMRES version of pseudo-arclength continuation [5], fixing the step 
in arclength to ds = .02, using a secant predictor [9], and beginning the continuation at 
c = 0, where the // = 1 is the solution. The vector with coordinates all equal to one is 
the solution of the discrete problem as well. We discretize the integral with the composite 
midpoint rule using 400 nodes. 

In Figure EISl we plot the average number of GMRES iterations per Newton iteration as 
a function of c. As one moves further on the path, the predictor becomes less effective, and 
the number of Newton iterations increase. Moreover, the norm of the solution also increases 
adding roughly one to the number of Krylov's per Newton. 
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